#################################################################################### 
####################################################################################
####################################################################################
# Appendix F & K: Figures F.3-F.5, K.6
####################################################################################
####################################################################################
library(here)
library(kableExtra)
library(lmtest)
library(multiwayvcov)
library(ggplot2)
####################################################################################
####################################################################################
## Read in survey data (from: analysis_results.R)
dat = readRDS(here("data/dat_unrestricted.rds"))
datb = readRDS(here("data/dat_restricted.rds"))

####################################################################################
####################################################################################
## Figures F.3

####################################################################################
# Health Care Groups Plot

d2 <- dat[!is.na(dat$health_care_groups),] %>% 
  group_by(treatment,health_care_groups) %>% 
  summarise(count=n()) %>% 
  mutate(meanc=mean(dat[!is.na(dat$health_care_groups) & dat$treatment == 0,]$health_care_groups)) %>%
  mutate(meant=mean(dat[!is.na(dat$health_care_groups) & dat$treatment == 1,]$health_care_groups)) %>%
  mutate(perc=count/sum(count)) %>%
  mutate(sample="Full (n=1477)")

d2b <- datb[!is.na(datb$health_care_groups),] %>% 
  group_by(treatment,health_care_groups) %>% 
  summarise(count=n()) %>% 
  mutate(meanc=mean(datb[!is.na(datb$health_care_groups) & datb$treatment == 0,]$health_care_groups)) %>%
  mutate(meant=mean(datb[!is.na(datb$health_care_groups) & datb$treatment == 1,]$health_care_groups)) %>%
  mutate(perc=count/sum(count)) %>%
  mutate(sample="Restricted (n=547)")

d2 = rbind(d2,d2b)

d2 = ggplot(d2, aes(x = factor(health_care_groups), y = perc*100, fill = factor(treatment))) +
  geom_bar(stat="identity", position = "dodge") +
  geom_text(aes(x=4, y=40, label=paste0("control:", "mu =","=", round(meanc, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  geom_text(aes(x=4, y=38, label=paste0("treatment:", "mu =","=", round(meant, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  theme_bw() + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(d2$perc*100)+1), labels=function(x) paste0(x,"%")) +
  scale_x_discrete(labels= c("(1) A\nStrongly", "(2) A", "(3)\nNeither",
                             "(4) B\n", "(5) B\nStrongly")) +
  labs(title = "Preferred Role of Government and NGOs in Health Care Provision",
       x = "(A) Govt provides most care (B) NGOs provide most care",
       y = NULL) + 
  scale_fill_manual(values=c("#999999", "#E69F00"),
                    labels=c("Control", "Treatment")) +
  theme(legend.title=element_blank(),
        plot.title = element_text(size=14,hjust = 0.5),
        #text=element_text(family="CM Roman"),
        legend.position = c(0.9, 0.9),
        legend.background = element_blank()) +
  facet_wrap(~sample)

ggsave(file = here("plots/expectation_group_wrap.pdf", sep = ""), d2, width = 8, height = 5)

rm(d2,d2b)

####################################################################################
####################################################################################
## Figures F.4

####################################################################################
# Health Care Payment Plot
d2 <- dat[!is.na(dat$health_care_pay),] %>% 
  group_by(treatment,health_care_pay) %>% 
  summarise(count=n()) %>% 
  mutate(meanc=mean(dat[!is.na(dat$health_care_pay) & dat$treatment == 0,]$health_care_pay)) %>%
  mutate(meant=mean(dat[!is.na(dat$health_care_pay) & dat$treatment == 1,]$health_care_pay)) %>%
  mutate(perc=count/sum(count)) %>%
  mutate(sample="Full (n=1477)")

d2b <- datb[!is.na(datb$health_care_pay),] %>% 
  group_by(treatment,health_care_pay) %>% 
  summarise(count=n()) %>% 
  mutate(meanc=mean(datb[!is.na(datb$health_care_pay) & datb$treatment == 0,]$health_care_pay)) %>%
  mutate(meant=mean(datb[!is.na(datb$health_care_pay) & datb$treatment == 1,]$health_care_pay)) %>%
  mutate(perc=count/sum(count)) %>%
  mutate(sample="Restricted (n=547)")

d2 = rbind(d2,d2b)

d2 = ggplot(d2, aes(x = factor(health_care_pay), y = perc*100, fill = factor(treatment))) +
  geom_bar(stat="identity", position = "dodge") +
  geom_text(aes(x=5.5, y=35, label=paste0("control:", "mu =","=", round(meanc, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  geom_text(aes(x=5.5, y=33, label=paste0("treatment:", "mu =","=", round(meant, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  theme_bw() + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(d2$perc*100)+1), labels=function(x) paste0(x,"%")) +
  scale_x_discrete(labels= c("(1) A\nStrongly", "(2) A", "(3) B", 
                             "(4) B\nStrongly", "(5) C", "(6) C\nStrongly")) + 
  labs(title = "Preferred Role of Government and NGOs in Health Care Financing and Provision",
       x = "(A) Govt provides and pays (B) Govt pays, NGOs provide (C) NGOs provide and pay",
       y = NULL) + 
  scale_fill_manual(values=c("#999999", "#E69F00"),
                    labels=c("Control", "Treatment")) +
  theme(legend.title=element_blank(),
        plot.title = element_text(size=14,hjust = 0.5),
        #text=element_text(family="CM Roman"),
        legend.position = c(0.9, 0.9),
        legend.background = element_blank()) +
  facet_wrap(~sample)

ggsave(file = here("plots/expectation_pay_wrap.pdf", sep = ""), d2, width = 8, height = 5)

rm(d2,d2b)

####################################################################################
####################################################################################
## Figures F.5

####################################################################################
# Donation Plot
d2 <- dat[!is.na(dat$donate_lg),] %>% 
  group_by(treatment,donate_lg) %>% 
  summarise(count=n()) %>% 
  mutate(meanc=mean(dat[!is.na(dat$donate_lg) & dat$treatment == 0,]$donate_lg)) %>%
  mutate(meant=mean(dat[!is.na(dat$donate_lg) & dat$treatment == 1,]$donate_lg)) %>%
  mutate(perc=count/sum(count)) %>%
  mutate(sample="Full (n=1477)")

d2b <- datb[!is.na(datb$donate_lg),] %>% 
  group_by(treatment,donate_lg) %>% 
  summarise(count=n()) %>% 
  mutate(meanc=mean(datb[!is.na(datb$donate_lg) & datb$treatment == 0,]$donate_lg)) %>%
  mutate(meant=mean(datb[!is.na(datb$donate_lg) & datb$treatment == 1,]$donate_lg)) %>%
  mutate(perc=count/sum(count)) %>%
  mutate(sample="Restricted (n=547)")

a = d2b[19:20,]
a$donate_lg = c(8,9)
a$perc = c(0,0)

d2 = rbind(d2,d2b,a)


donation = ggplot(d2, aes(x = factor(donate_lg), y = perc*100, fill = factor(treatment))) +
  geom_bar(stat="identity", position = "dodge") +
  theme_bw() + 
  geom_text(aes(x=9, y=25, label=paste0("control:", "mu =","=", round(meanc, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  geom_text(aes(x=9, y=23, label=paste0("treatment:", "mu =","=", round(meant, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(d2$perc*100)+1), labels=function(x) paste0(x,"%")) +
  scale_x_discrete(labels= c("0\nAll Coins\nto Govt","1","2","3","4","5","6","7","8","9","10\nAll Coins\nto NGO")) + 
  labs(title = "Distribution of Donation Allocation between Government and CHP Program",
       x = "Number of Coins Donated to NGO",
       y = NULL) + 
  scale_fill_manual(values=c("#999999", "#E69F00"),
                    labels=c("Control", "Treatment")) +
  theme(legend.title=element_blank(),
        plot.title = element_text(size=14,hjust = 0.5),
        #text=element_text(family="CM Roman"),
        legend.position = c(0.9, 0.9),
        legend.background = element_blank()) +
  facet_wrap(~sample)

ggsave(file = here("plots/expectation_donation_wrap.pdf", sep = ""), donation, width = 8, height = 5)

rm(d2,d2b)

####################################################################################
####################################################################################
## Figures K.6

####################################################################################
# Policy Priorities
d2 <- dat[!is.na(dat$priority_health_nat),] %>% 
  group_by(treatment,priority_health_nat) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count)) %>%
  mutate(meanc=mean(dat[!is.na(dat$priority_health_nat) & dat$treatment == 0,]$priority_health_nat)) %>%
  mutate(meant=mean(dat[!is.na(dat$priority_health_nat) & dat$treatment == 1,]$priority_health_nat)) %>%
  mutate(sample="Full (n=1477)")


d2b <- datb[!is.na(datb$priority_health_nat),] %>% 
  group_by(treatment,priority_health_nat) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count)) %>%
  mutate(meanc=mean(datb[!is.na(datb$priority_health_nat) & datb$treatment == 0,]$priority_health_nat)) %>%
  mutate(meant=mean(datb[!is.na(datb$priority_health_nat) & datb$treatment == 1,]$priority_health_nat)) %>%
  mutate(sample="Restricted (n=547)")

d2 = rbind(d2,d2b)

health_priority = ggplot(d2, aes(x = factor(priority_health_nat), y = perc*100, fill = factor(treatment))) +
  geom_bar(stat="identity", position = "dodge") +
  theme_bw() + 
  geom_text(aes(x=2, y=30, label=paste0("control:", "mu =","=", round(meanc, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  geom_text(aes(x=2, y=28, label=paste0("treatment:", "mu =","=", round(meant, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(d2$perc*100)+1), labels=function(x) paste0(x,"%")) +
  scale_x_discrete(labels= c("Not\nListed","Third\nPriority","Second\nPriority","First\nPriority")) +
  labs(title = "Prioritization of Health Problem for Government to Address",
       x = "Prioritization of Health",
       y = NULL) + 
  scale_fill_manual(values=c("#999999", "#E69F00"),
                    labels=c("Control", "Treatment")) +
  theme(legend.title=element_blank(),
        plot.title = element_text(size=14,hjust = 0.5),
        #text=element_text(family="CM Roman"),
        legend.position = c(0.1, 0.9),
        legend.background = element_blank()) +
  facet_wrap(~sample)

ggsave(file = here("plots/health_priority_wrap.pdf", sep = ""), health_priority, width = 8, height = 5)
####################################################################################